NSF-AST and DEB topic modeling - head

Useful tools… to find NA in a column of a data frame: new_DF <- DF[is.na(DF$Var),] eg. new_lemma_new <- lemma_new[is.na(lemma_new$word),]

library(styler)

#### Functions #####

#--------------------
readpercent <- function()
{ 
  n <- readline(prompt="Enter a percent, positive number selects head & negative selects tail[-100-100]: ")
  if(!(grepl("^-?[0-9]+$",n)) || as.integer(n) > 100 || as.integer(n) < -100) #
  {
    return(readpercent())
  }
  return(as.integer(n))
}

#--------------------
HistogramFun<-function(funding, fstatedPercent, numRows) {
  cat("typeof funding",typeof(funding), "typeof nrows", typeof(numRows), "nrows funding ", typeof(numRows), "numRows::", numRows, "\n")
  if (fstatedPercent < 0) {
    color1 <- "gray"; 
    color2 <- "gray40";
    nInHead <- as.integer(((100 + fstatedPercent)/100) * numRows)
    nInTail <- numRows - nInHead
    headPercent <- fstatedPercent + 100
    tailPercent <- abs(fstatedPercent)
  } else {
    color2 <- "gray"
    color1 <- "gray40"
    nInHead <- as.integer((fstatedPercent/100) * numRows)
    nInTail <- numRows - nInHead
    headPercent <- fstatedPercent
    tailPercent <- 100 - fstatedPercent
  }
  cols = c(rep(color1, times = nInHead), rep(color2, times = nInTail))
  #ylim= c(0, 3),  
  barplot(funding, col = cols, border=cols, main = NULL, 
          xpd = FALSE, las=1
  )
  
  #legend("right", legend = c("Head (10%)", "Long Tail (90%)"), fill = c("gray40", "gray"), cex=.9)
  legend("right", 
         legend = c(paste("Head (", headPercent, "%)", sep = ""), 
                    paste("Tail (", tailPercent, "%)", sep = "")), 
         fill = c(color1, color2), cex=.9)
  title("The focus for later analysis is shaded dark", cex.main=.95, line= -2) 
  title(xlab="Size Ordered Grants - bar for each grant", line=0)
  title(ylab="Dollar Value of Award in Millions", line = 2) 
}

# Read in the file
MyData <- read.csv("data/CNS.csv", header=TRUE, encoding = "UTF-8", stringsAsFactors=FALSE) # tab delimited

# pull off the directory part of the file name to use later for saving output files.
#topicDir <- tools::file_path_sans_ext(filename)
topicDir <- "Awards_CNS"

dir.create(topicDir, showWarnings = FALSE)



# decide on which parts of the curve to analyze
# User must enter the # on the console!
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

nsf_funding <- data_frame(id = as.character(MyData$AwardNumber),
                          abstract = MyData$Abstract,
                          amount = MyData$AwardedAmountToDate)

Warning: `data_frame()` was deprecated in tibble 1.1.0.
Please use `tibble()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

# remove html tags in the form <anything>
# ??? drop records with empty abstracts
# changed PBH 2019-02-27 to put in a space when html is replaced
nsf_funding$abstract <- gsub("<[^>]+>", " ", nsf_funding$abstract)
# remove where the word is a number
nsf_funding$abstract <- gsub("\\d", "", nsf_funding$abstract)
# if there are any "words" with imbedded periods Eg "John.P.Doud, replace the peroid with a space.
nsf_funding$abstract <- gsub("\\.", " ", nsf_funding$abstract)
# remove "'s"
nsf_funding$abstract <- gsub("'s", "", nsf_funding$abstract) # to fix eg trichoptera:hydropsychidae 
nsf_funding$abstract <- gsub(":", " ", nsf_funding$abstract)
#nsf_funding$abstract <- nsf_funding[-grep('^\\d+$', nsf_funding$abstract),]
# Convert the string repreentation of numbers in the amount field by removing "$" and ","
nsf_funding$amount <- as.numeric(gsub('[$,]', '', nsf_funding$amount))

# sort the abstracts by the amount awarded
#abstract_sort <- nsf_funding %>% 
nsf_funding <- nsf_funding[order(-nsf_funding$amount), ]

statedPercent <- 100
# if statedPercent is not set in the code ask the user
if(!exists("statedPercent")){statedPercent <- readpercent();}
# if it is negative then it is really the tail and the head is the length the array minus that amount

#This displays on the screen
HistogramFun(nsf_funding$amount/1000000, statedPercent, nrow(nsf_funding))

typeof funding double typeof nrows integer nrows funding  integer numRows:: 3000 

#This saves a file
png(paste(topicDir,"/focusHistogram.png", sep=""), width = 5, height = 4, units = 'in', res = 300)

HistogramFun(nsf_funding$amount/1000000, statedPercent, nrow(nsf_funding))

typeof funding double typeof nrows integer nrows funding  integer numRows:: 3000 

dev.off()

quartz_off_screen 
                2 

if (statedPercent != 100) {
  # Only do this if you want to only analyse part of the collection
  nsf_funding <- head(nsf_funding, as.integer(nrow(nsf_funding)*(headpercent/100)))
}

For astronomy data may need to hide the few big projects to see the curve.

### Only run this to use custom stop words. ####
library(tidytext)
#Custom stop words
#adapt for Bio!!!
my_stop_words <- bind_rows(stop_words, 
                           data_frame(word = c("formation", "science", "stars", "galaxies",
                                               "data", "star",
                                               "project", "astronomy", "research", "students", "galaxy",
                                               "program", "public", "observations", "study",
                                               "graduate", "results", "team", "award",
                                               "dr", as.character(1:12)),
                                                
                                      lexicon = rep("custom", 32)))

library(tidyverse) # R is better when it is tidy

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.3     ✓ purrr   0.3.4
✓ tibble  3.0.4     ✓ stringr 1.4.0
✓ tidyr   1.1.2     ✓ forcats 0.5.0
✓ readr   1.4.0     
── Conflicts ────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(stringr)  # for string manipulation
# Save copy of abstracts for bigram work
library(tidytext)
nsf_funding_bigram <- nsf_funding
nsf_funding_save <- nsf_funding # This will be used to find bigrams in the documents.

nsf_funding <- nsf_funding %>% 
  unnest_tokens(word, abstract)  %>% 
  #anti_join(my_stop_words) # Add custom stop words.
  anti_join(stop_words) # Use normal stop words.

Joining, by = "word"

# remove rows that are null.
#nsf_funding <- nsf_funding[!(nsf_funding$word == null), ]
nsf_funding <- nsf_funding[!is.null(nsf_funding$word) ]
nsf_funding %>% filter(nsf_funding$word != "null" )

nsf_funding_bigram <- nsf_funding_bigram %>% 
  unnest_tokens(bigram, abstract, token = "ngrams", n = 2) %>% 
  #anti_join(my_stop_words) # Add custom stop words.
  anti_join(stop_words, by = c("bigram" = "word")) # Use normal stop words.
lemma_unique <- nsf_funding %>%
  select(word) %>%
  mutate(word_clean = str_replace_all(word,"\u2019s|'s","")) %>%
#  mutate(word_clean = ifelse(str_detect(word_clean,"[^[:alpha:]]"),NA,word_clean)) %>%
  # this :alpha: causes problems with UTF-8
 # mutate(word_clean = ifelse(str_detect(word_clean,"[^[:alpha:]]"),NA,word_clean)) %>%
  filter(!duplicated(word_clean)) %>%
  filter(!is.na(word_clean)) %>%
  arrange(word)

#load packages koRupus and sylly
library(dplyr)
library(sylly)
library(koRpus)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
For information on available language packages for 'koRpus', run

  available.koRpus.lang()

and see ?install.koRpus.lang()


Attaching package: ‘koRpus’

The following object is masked from ‘package:readr’:

    tokenize

#install.koRpus.lang("en")
library(koRpus.lang.en)
# setup for lemmanization
# Outside of R Install tree-tagger from https://cran.r-project.org/web/packages/koRpus/vignettes/koRpus_vignette.html
# http://www.cis.uni-muenchen.de/~schmid/tools/TreeTagger/#parfiles
# to install tree-tagger in order to run koRpus in R
system.time(
  lemma_tagged <- treetag(lemma_unique$word_clean, treetagger="manual", 
                          
                          format="obj", TT.tknz=FALSE , lang="en",
                          TT.options=list(
                          path="Applications/TreeTagger", preset="en"),
                           debug=TRUE
                          )
)

        split=[[:space:]]
        ign.comp=-
        heuristics=abbr
        heur.fix=c("’", "'"), c("’", "'")
        sentc.end=., !, ?, ;, :
        detect=FALSE, FALSE
        clean.raw=
        perl=FALSE
        stopwords=
        stemmer=
 
          TT.tokenizer:  koRpus::tokenize() 
          tempfile: /var/folders/gp/rs2l_qcn4nld9bwwl_c5y5wh0000gn/T//Rtmpq07u1C/tokenize2056224d1f7b.txt 
          file:  /var/folders/gp/rs2l_qcn4nld9bwwl_c5y5wh0000gn/T//Rtmpq07u1C/tempTextFromObject20565d2ea8d0.txt 
          TT.lookup.command:   
          TT.pre.tagger:  grep -v '^$' | 
          TT.tagger:  Applications/TreeTagger/bin/tree-tagger 
          TT.opts:  -token -lemma -sgml -pt-with-lemma -quiet 
          TT.extra.opts:   
          TT.params:  /Users/zhengsongyi/Desktop/2022Spring/Capstone/test/Applications/TreeTagger/lib/english.par 
          TT.filter.command:  | perl -pe 's/\tV[BDHV]/\tVB/;s/IN\/that/\tIN/;' 

          sys.tt.call:  cat  /var/folders/gp/rs2l_qcn4nld9bwwl_c5y5wh0000gn/T//Rtmpq07u1C/tokenize2056224d1f7b.txt |  grep -v '^$' | Applications/TreeTagger/bin/tree-tagger -token -lemma -sgml -pt-with-lemma -quiet  /Users/zhengsongyi/Desktop/2022Spring/Capstone/test/Applications/TreeTagger/lib/english.par | perl -pe 's/\tV[BDHV]/\tVB/;s/IN\/that/\tIN/;' 

Assuming 'UTF-8' as encoding for the input file. If the results turn out to be erroneous, check the file for invalid characters, e.g. em.dashes or fancy quotes, and/or consider setting 'encoding' manually.

   user  system elapsed 
  0.998   0.087   1.049 

# Make sure all packages including tidyverse is loaded or you may get an error that the preset is not recognized. In particular check koRpus.lang.en

#lemma_tagged_tbl <- tbl_df(lemma_tagged@TT.res)
#lemma_tagged_tbl <- tibble(lemma_tagged@TT.res)
lemma_tagged_tbl <- as_tibble(lemma_tagged)
# replace unknowns with the old raw word
#lemma_unique <- lemma_unique %>% 
#  left_join(lemma_tagged_tbl %>%
#              filter(lemma != "<unknown>") %>%
#              select(token, lemma, wclass), by = c("word_clean" = "token")
#            ) %>%
#  arrange(word)
# replace wclass w/ tag
lemma_unique <- lemma_unique %>% 
  left_join(lemma_tagged_tbl %>%
              filter(lemma != "<unknown>") %>%
              select(token, lemma, tag), by = c("word_clean" = "token")
            ) %>%
  arrange(word)

# Replace all NA with the word in word_clean
lemma_unique = within(lemma_unique, {
  lemma = ifelse(is.na(lemma), word_clean, lemma)  } )

#replace all of the original words in the nsf_funding with lemmas
# This also needs to be done for the bigrams where the substatution happens for each of the two columns in the bigram. ???
lemma_new <- lemma_unique
#lemma_new <- left_join(nsf_funding,lemma_new, by = "word")
lemma_new <- left_join(nsf_funding,lemma_new, by = "word") %>%
  mutate(word = ifelse(is.na(lemma), word_clean, lemma)) # %>%
lemma_new <- subset(lemma_new, select = c(id, word))
lemma_new[is.na(lemma_new$word),]

NA

There should be no NA words!!! if so there is a problem with the input… UTF conversion problem.

library(tm)

Loading required package: NLP

Attaching package: ‘NLP’

The following object is masked from ‘package:ggplot2’:

    annotate


Attaching package: ‘tm’

The following object is masked from ‘package:koRpus’:

    readTagged

#myCorpus <- Corpus(VectorSource(nsf_funding$word)) 
myCorpus <- Corpus(VectorSource(lemma_new$word)) 

#create the toSpace content transformer
toSpace <- content_transformer(function(x, pattern) {return (gsub(pattern, " ", x))})

#Strip digits (std transformation, so no need for content_transformer)
myCorpus <- tm_map(myCorpus, removeNumbers)

Warning in tm_map.SimpleCorpus(myCorpus, removeNumbers) :
  transformation drops documents

dtm <- DocumentTermMatrix(myCorpus)

freq <- colSums(as.matrix(dtm))

Error: vector memory exhausted (limit reached?)

# Plot the occurence frequencies
wf=data.frame(term=names(freq),occurrences=freq)
library(ggplot2)
p <- ggplot(subset(wf, freq>20), aes(term, occurrences))
p <- p + geom_bar(stat="identity")
p <- p + theme(axis.text.x=element_text(angle=45, hjust=1))
p

#wordcloud
library(wordcloud)

Loading required package: RColorBrewer

#setting the same seed each time ensures consistent look across clouds
set.seed(42)
png(paste(topicDir,"/FreqWordCloud.png", sep=""), width=1280,height=800)
#limit words by specifying min frequency
wordcloud(names(freq),freq, max.words=100, scale=c(4,.2), colors=brewer.pal(6,"Dark2"))
while (!is.null(dev.list()))  dev.off()
#Second call puts the cloud image on the screen
wordcloud(names(freq), freq, max.words=100, scale=c(4,.2), colors=brewer.pal(6,"Dark2"))

Warning in wordcloud(names(freq), freq, max.words = 100, scale = c(4, 0.2),  :
  project could not be fit on page. It will not be plotted.

##Topic Modeling using lemmas

#word_counts <- nsf_funding %>%
word_counts <- lemma_new %>%
  #anti_join(my_stop_words) %>%
  count(id, word, sort = TRUE) %>%
  ungroup()

word_counts

#Remove numbers
#word_counts<-word_counts[-grep("\\b\\d+\\b", word_counts$word),]

desc_dtm <- word_counts %>%
  cast_dtm(id, word, n)

desc_dtm

<<DocumentTermMatrix (documents: 38, terms: 1752)>>
Non-/sparse entries: 4743/61833
Sparsity           : 93%
Maximal term length: 17
Weighting          : term frequency (tf)

#Calculate tf_idf
#word_counts <- word_counts %>%
 #bind_tf_idf(id, word, n)
#word_counts

#Continue tf_idf
#word_counts %>%
  #arrange(tf_idf)

#Calculate tf-idf
#desc_tf_idf <- nsf_funding %>% 
desc_tf_idf <- lemma_new %>% 
  count(id, word, sort = TRUE) %>%
  ungroup() %>%
  bind_tf_idf(word, id, n)

#New document term matrix for topic modeling with custom stop words #??? rather than low idf use low tf_idf - may be replaced with a trim of the word list by ft-idf low_idf <- subset(desc_tf_idf, idf>0.38 & idf<2) #topicFileName <- gsub(“.csv”,“low_idf.csv”, filename) write.csv(low_idf,paste(topicDir,“/low_idf.csv”, sep=““))

#What are highest tf-idf words?
desc_tf_idf %>% 
  arrange(-tf_idf)

#What are lowest idf words?
desc_tf_idf %>% 
  arrange(idf)

freq <- colSums(as.matrix(desc_dtm))

#create sort order (descending)
ord <- order(freq,decreasing=TRUE)

#inspect most frequently occurring terms
freq[head(ord)]

  project     datum  research    system algorithm   network 
      116       101        82        80        74        66 

#inspect least frequently occurring terms
freq[tail(ord)]

transaction     tunable      update      usable        weak  workaround 
          1           1           1           1           1           1 

#Ignore this maybe??? do not use.
#wf=data.frame(term=names(freq),occurrences=freq)
#library(ggplot2)
#topicFileName <- gsub(".csv","tf-idfhistogram.png", filename)
#p <- ggplot(subset(wf, freq>250), aes(term, occurrences))
#p <- p + geom_bar(stat="identity")
#p <- p + theme(axis.text.x=element_text(angle=45, hjust=1))
#p
#ggsave(topicFileName, device = "png")

#wordcloud
library(wordcloud)
#setting the same seed each time ensures consistent look across clouds
set.seed(42)
png(paste(topicDir,"/TF-IDFWordCloud.png", sep=""), width=1280,height=800)
#limit words by specifying min frequency
#wordcloud(names(freq),freq, min.freq=80,colors=brewer.pal(6,"Dark2"))
wordcloud(names(freq),freq, min.freq=20,scale=c(4,.2), colors=brewer.pal(6,"Dark2"))
while (!is.null(dev.list()))  dev.off()
wordcloud(names(freq),freq, min.freq=20,scale=c(4,.2), colors=brewer.pal(6,"Dark2"))

Warning in wordcloud(names(freq), freq, min.freq = 20, scale = c(4, 0.2),  :
  project could not be fit on page. It will not be plotted.

https://cran.r-project.org/web/packages/ldatuning/vignettes/topics.html Select number of topics for LDA model Murzintcev Nikita 2016-10-24

Package can be installed from CRAN Package ldatuning realizes 4 metrics to select perfect number of topics for LDA model.

#LDA Tuning: Pick the Number of topics

Warning message:
restarting interrupted promise evaluation 

#install.packages("ldatuning")
# Changed method in the exmaple from "Gibbs"" to "VEM" to mathc the Nasa Example
library("ldatuning")
library(topicmodels)
system.time(
  result <- FindTopicsNumber(
  desc_dtm,
  topics = seq(from = 10, to = 50, by = 1),
  metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),
  method = "Gibbs",
  control = list(seed = 1234),
  mc.cores = 4L,
  verbose = TRUE
))

fit models... done.
calculate metrics:
  Griffiths2004... done.
  CaoJuan2009... done.
  Arun2010... done.
  Deveaud2014... done.
   user  system elapsed 
  3.240   0.367  41.289 

png(paste(topicDir,"/FindTopicNumber.png", sep=""), width=1280,height=800)
FindTopicsNumber_plot(result)

while (!is.null(dev.list()))  dev.off()
FindTopicsNumber_plot(result)

Gibbs method consensus shows a peak of 23 for the DEB data Gibbs method consensus shows a peak of 30 for the AST data

harmonicMean <- function(logLikelihoods, precision = 2000L) {
  llMed <- median(logLikelihoods)
  as.double(llMed - log(mean(exp(-mpfr(logLikelihoods,
                                       prec = precision) + llMed))))
}

seqk <- seq(2, 100, 1)
burnin <- 1000
iter <- 1000
keep <- 50
system.time(fitted_many <- lapply(seqk, function(k) topicmodels::LDA(desc_dtm, k = k,
                                                     method = "Gibbs",control = list(burnin = burnin,
                                                                         iter = iter, keep = keep) )))

   user  system elapsed 
299.010   1.881 302.937 

roundto10 <- function(x) ceiling(max(x)/10)*10

# extract logliks from each topic
library("Rmpfr")

Loading required package: gmp

Attaching package: ‘gmp’

The following objects are masked from ‘package:base’:

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: ‘Rmpfr’

The following object is masked from ‘package:gmp’:

    outer

The following objects are masked from ‘package:stats’:

    dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm

The following objects are masked from ‘package:base’:

    cbind, pmax, pmin, rbind

logLiks_many <- lapply(fitted_many, function(L)  L@logLiks[-c(1:(burnin/keep))])

# compute harmonic means
hm_many <- sapply(logLiks_many, function(h) harmonicMean(h))

#Now that is have calculated the harmonic means of the models for 2 - 100 topics, you can plot the results using ggplot2 package.

require(scales)

Loading required package: scales

Attaching package: ‘scales’

The following object is masked from ‘package:purrr’:

    discard

The following object is masked from ‘package:readr’:

    col_factor

ldaplot <- ggplot(data.frame(seqk, hm_many), aes(x=seqk, y=hm_many)) + geom_path(lwd=1.5) +
  theme(text = element_text(family= NULL),
        axis.title.y=element_text(vjust=1, size=16),
        axis.title.x=element_text(vjust=-.5, size=16),
        axis.text=element_text(size=16),
        plot.title=element_text(size=20)) +
  scale_y_continuous(labels = comma) +
  xlab('Number of Topics') +
  ylab('Harmonic Mean') +
  ggplot2::annotate("text", x = 25, y = roundto10(max(hm_many)+10000), # calculate the max rather than hard code ???
           label = paste("The optimal number of topics is", seqk[which.max(hm_many)])) +
  ggtitle(expression(atop("LDA Analysis of NSF Program", 
                          atop(italic("Harmonic Mean"), ""))))
 
png(paste(topicDir,"/LDAHamonicMean.png", sep=""), width=1280,height=800)
ldaplot
while (!is.null(dev.list()))  dev.off()
ldaplot

#Harmonic mean above has optimum of 20 for combined collaborations in astronomy. Prior analysis used 22 but the peak is about the same.
#Inspired by: https://www.tidytextmining.com/nasa.html#topic-modeling
library(topicmodels)

#Harmonic mean was 34 topics but 29 is almost the same; 
# [23] -408680.1 -407910.1 -408071.9 -407422.7 -407173.9 -406342.0 -407583.3 -407614.1 -407782.2 -408115.3 -405507.6
# be aware that running this model is time intensive
NumTopics <- 18
desc_lda <- LDA(desc_dtm, method = "Gibbs", k = NumTopics, control = list(seed = 1234))
desc_lda

A LDA_Gibbs topic model with 18 topics.

# Since it is time consuming to create the desc_lda file we save it here for backup. 
saveRDS(desc_lda, file = paste(topicDir,"/desc_lda.rds", sep = ""))
# to get it back desc_lda <- load(paste(topicDir,"/desc_lda.RData"))

library(tidytext)
tidy_lda <- tidy(desc_lda)

tidy_lda

#Examine top 10 terms for each topic
top_terms <- tidy_lda %>%
  group_by(topic) %>%
  top_n(10, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)

top_terms

#https://tm4ss.github.io/docs/Tutorial_6_Topic_Models.html

top10termsPerTopic <- terms(desc_lda, 10)
topicNames <- apply(top10termsPerTopic, 2, paste, collapse=" ")

# have a look at some of the results (posterior distributions)
tmResult <- posterior(desc_lda)

library(wordcloud)
# visualize topics as word cloud
#for (TopicNum in 1:NumTopics){
for (TopicNum in 1:18){
  set.seed(42)
  topicToViz <- TopicNum # change for your own topic of interest
  #topicToViz <- 1 # change for your own topic of interest
  #topicToViz <- grep('telescope', topicNames)[1] # Or select a topic by a term contained in its name
  # select to 100 most probable terms from the topic by sorting the term-topic-probability vector in decreasing order
  top100terms <-
  sort(tmResult$terms[topicToViz, ], decreasing = TRUE)[1:100]
  words <- names(top100terms)
  # extract the probabilites of each of the 100 terms
  probabilities <-
  sort(tmResult$terms[topicToViz, ], decreasing = TRUE)[1:100]
  # visualize the terms as wordcloud
  mycolors <- brewer.pal(8, "Dark2")
  #topicFileName <- gsub(".csv", paste("Topic", toString(TopicNum), "WordCloud.png", sep = ""), filename) 
  png(paste(topicDir,"/WordCloud.png", sep=""), width = 1280, height = 800)
  while (!is.null(dev.list()))  dev.off()
  wordcloud(
  words,
  probabilities,
  scale = c(4, .2),
  random.order = FALSE,
  color = mycolors
  )
  dev.off()
  wordcloud(
  words,
  probabilities,
  scale = c(4, .2),
  random.order = FALSE,
  color = mycolors
  )
}

Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  development could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  accurately could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  computationally could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  discover could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  primitive could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  modeling could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  difficult could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  source could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  compact could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  threshold could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  combinatorial could not be fit on page. It will not be plotted.

Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  computational could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  accurately could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  computationally could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  transcriptome could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  speed could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  proposal could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  compact could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  threshold could not be fit on page. It will not be plotted.
Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  optimal could not be fit on page. It will not be plotted.

Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  communication could not be fit on page. It will not be plotted.

Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  communication could not be fit on page. It will not be plotted.

Warning in wordcloud(words, probabilities, scale = c(4, 0.2), random.order = FALSE,  :
  performance could not be fit on page. It will not be plotted.

library(ggplot2)
png(paste(topicDir,"/TermsLDAHistogram.png", sep=""), width = 1280, height = 1600)
top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  group_by(topic, term) %>%    
  arrange(desc(beta)) %>%  
  ungroup() %>%
  mutate(term = factor(paste(term, topic, sep = "__"), 
                       levels = rev(paste(term, topic, sep = "__")))) %>%
  ggplot(aes(term, beta, fill = as.factor(topic))) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  scale_x_discrete(labels = function(x) gsub("__.+$", "", x)) +
  labs(title = "Top terms in each LDA topic",
       x = NULL, y = expression(beta)) +
  facet_wrap(~ topic, ncol = 3, scales = "free")
#topicFileName <- gsub(".csv","TermsLDAHistogram.png", filename)
#ggsave(topicFileName, height = 14, device = "png")
#png(paste(topicDir,"/TermsLDAHistogram.png", sep=""), width = 1280, height = 1600)
while (!is.null(dev.list()))  dev.off()

# Calculate the top terms in each topic
top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  group_by(topic, term) %>%    
  arrange(desc(beta)) %>%  
  ungroup() %>%
  mutate(term = factor(paste(term, topic, sep = "__"), 
                       levels = rev(paste(term, topic, sep = "__")))) %>%
  ggplot(aes(term, beta, fill = as.factor(topic))) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  scale_x_discrete(labels = function(x) gsub("__.+$", "", x)) +
  labs(title = "Top terms in each LDA topic",
       x = NULL, y = expression(beta)) +
  facet_wrap(~ topic, ncol = 3, scales = "free")

#Top 10 terms in each model

terms(desc_lda, 10)

      Topic 1      Topic 2      Topic 3      Topic 4       Topic 5        Topic 6      Topic 7     Topic 8    Topic 9       
 [1,] "method"     "datum"      "learn"      "student"     "linear"       "design"     "project"   "design"   "compute"     
 [2,] "develop"    "cluster"    "training"   "software"    "algebraic"    "auction"    "analysis"  "base"     "application" 
 [3,] "approach"   "privacy"    "research"   "research"    "method"       "resource"   "technique" "software" "parallel"    
 [4,] "datum"      "machine"    "algorithm"  "reliability" "project"      "guarantee"  "tool"      "system"   "architecture"
 [5,] "gene"       "address"    "run"        "provide"     "construction" "allocation" "market"    "approach" "emerge"      
 [6,] "tool"       "analysis"   "connection" "science"     "tool"         "mechanism"  "algorithm" "analysis" "program"     
 [7,] "result"     "future"     "efficient"  "university"  "dimension"    "objective"  "field"     "tradeoff" "research"    
 [8,] "expression" "discovery"  "task"       "include"     "pseudorandom" "require"    "question"  "space"    "fault"       
 [9,] "cell"       "framework"  "improve"    "safety"      "polynomial"   "art"        "develop"   "level"    "gpu"         
[10,] "technology" "researcher" "machine"    "reu"         "set"          "agent"      "model"     "complex"  "include"     
      Topic 10      Topic 11        Topic 12      Topic 13        Topic 14    Topic 15        Topic 16      Topic 17     
 [1,] "network"     "algorithm"     "graph"       "theory"        "support"   "datum"         "system"      "system"     
 [2,] "develop"     "optimization"  "network"     "information"   "impact"    "network"       "project"     "efficiency" 
 [3,] "clock"       "focus"         "processing"  "study"         "time"      "communication" "device"      "energy"     
 [4,] "topology"    "framework"     "datum"       "project"       "broad"     "code"          "batteryless" "memory"     
 [5,] "performance" "student"       "social"      "science"       "signal"    "scheme"        "computer"    "performance"
 [6,] "project"     "submodular"    "challenge"   "understanding" "deep"      "research"      "outreach"    "propose"    
 [7,] "scientific"  "algorithmic"   "research"    "interactive"   "domain"    "caching"       "reliable"    "hardware"   
 [8,] "tree"        "modern"        "practice"    "community"     "structure" "practical"     "energy"      "scale"      
 [9,] "circuit"     "aim"           "cost"        "measurement"   "property"  "information"   "student"     "solution"   
[10,] "mode"        "undergraduate" "interaction" "compression"   "model"     "storage"       "industry"    "computer"   
      Topic 18        
 [1,] "language"      
 [2,] "programme"     
 [3,] "client"        
 [4,] "implementation"
 [5,] "component"     
 [6,] "interface"     
 [7,] "programmer"    
 [8,] "synthesis"     
 [9,] "program"       
[10,] "exist"         

##Examine which topics are associated with which documents

lda_gamma <- tidy(desc_lda, matrix = "gamma")

lda_gamma

library(data.table)

data.table 1.13.6 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

Attaching package: ‘data.table’

The following object is masked from ‘package:purrr’:

    transpose

The following objects are masked from ‘package:dplyr’:

    between, first, last

# Change the index to be a column. nsf_funding_save is in funding order.
setDT(nsf_funding_save, keep.rownames = TRUE)[]

#Change the column name form the "rn" defaculty to "FundingRank"
colnames(nsf_funding_save)[1] <- "FundingRank"
# Change class of nsf_funding_save$id to charcater(from integer)
nsf_funding_save$id <- as.character(nsf_funding_save$id)
lda_gamma2 <- dplyr::left_join(lda_gamma, nsf_funding_save, 
              select(FundingRank), by = c("document" = "id")) 
# convert FundingRank to a number
lda_gamma2$FundingRank = as.numeric(lda_gamma2$FundingRank)#[lda_gamma2]
# convert Topic to a character

# Sort the topics by the amont of monney they braught in. Sum of document Gamma = 1
# The funding associated with each topic is the 
# sum of all (gamma*funding amount for each grant)
lda_gamma2 <- lda_gamma2 %>% mutate(worth = gamma * amount)
# calculate the sum of the worth for each topic
topic_worth <- aggregate(lda_gamma2$worth, by=list(topic = lda_gamma2$topic), FUN = "sum")
#topic_worth <- aggregate(lda_gamma2$worth, by=list(topic = as.character(lda_gamma2$topic)), FUN = "sum")
colnames(topic_worth)[2] <- "WorthSum"
# sort by most worthy of topic
topic_worth <- topic_worth %>% arrange(desc(WorthSum))
# save topic_worth as CVS

setDT(topic_worth, keep.rownames = "TopicWorthRank")[]

lda_gamma3 <- left_join(lda_gamma2, topic_worth, 
              select(x, TopicWorthRank), by = "topic") 
lda_gamma3$TopicWorthRank = as.numeric(lda_gamma3$TopicWorthRank)

p <- ggplot(lda_gamma3, aes(x=FundingRank, y=TopicWorthRank)) +
    geom_point(aes(size=gamma), shape = 17) +   # draw points
    scale_size_continuous(limits=c(.2,1)) +
    xlim(c(0, length(nsf_funding_save$id))) + 
   ylim(c(1, NumTopics)) +
   labs(subtitle="Topic by Document", 
       y="Topic #", 
       x="Document #", 
       title="Gamma Level sorted by Funding Level", 
       caption="Topic Analysis")
p

Warning: Removed 648 rows containing missing values (geom_point).

ggsave(paste(topicDir,"/DocTopic.png", sep=""), height = 10, device = "png")

Saving 20 x 10 in image
Warning: Removed 648 rows containing missing values (geom_point).

#png(paste(topicDir,"/TermsLDAHistogram.png", sep=""), width = 1280, height = 1600)
while (!is.null(dev.list()))  dev.off()
p

Warning: Removed 648 rows containing missing values (geom_point).

#Join funding amount to data
nsf_funding_join <- nsf_funding_save
nsf_funding_join$document <- as.character(nsf_funding_join$id)

funding_gamma <- merge(lda_gamma, nsf_funding_join, by="document", all.x = TRUE)
  funding_gamma <- subset(funding_gamma, select=-c(id, abstract))
  
funding_gamma

NA

#Need to join lda_gamma with funding amount vector, that would be informative
doc_top <- funding_gamma %>%
  group_by(document) %>%
#  top_n(1, gamma) %>%
  top_n(5, gamma) %>%
  ungroup() %>%
#  arrange(gamma)
  arrange(document, desc(gamma))
write.csv(doc_top, paste(topicDir,"/topTopics5Documents.csv", sep=""))

doc_top

#Consensus document for each topic
document_topics <- doc_top %>%
  count(document, topic, amount) %>%
  group_by(topic) %>%
  top_n(1, n) %>%
  ungroup() %>%
  transmute(consensus = document, topic, amount)

document_topics

#See which documents misidentified, this is probably not useful for us, since documents can easily be in multiple topics
doc_top %>%
  inner_join(document_topics, by = "topic") %>%
  count(document, consensus)

p

Warning: Removed 684 rows containing missing values (geom_text).

NO NOT USE THIS CHUNK. Use the next one to make lda_json for use in LDAvis package

#In this chunk I was trying to use LDAvis package, but was unable to convert topicmodels output

#install.packages(“LDAvis”) #devtools::install_github(“cpsievert/LDAvisData”)

#’ Convert the output of a topicmodels Latent Dirichlet Allocation to JSON #’ for use with LDAvis #’ #’ @param fitted Output from a topicmodels model. = desc_lda #’ @param corpus Corpus object used to create the document term #’ matrix for the model. This should have been create with #’ the tm package’s function. = myCorpus #’ @param doc_term The document term matrix used in the #’ model. This should have been created with the tm package’s #’ function. = dtm #’ #’ @seealso . #’ @export

topicmodels_json_ldavis <- function(fitted, corpus, doc_term){ # Required packages library(topicmodels) library(dplyr) library(stringi) library(tm) library(LDAvis)

# Find required quantities

phi <- posterior(fitted)\(terms %>% as.matrix theta <- posterior(fitted)\)topics %>% as.matrix vocab <- colnames(phi) doc_length <- vector() for (i in 1:length(corpus)) { temp <- paste(corpus[[i]]$content, collapse = ’ ‘) doc_length <- c(doc_length, stri_count(temp, regex =’\S+’)) }

temp_frequency <- as.matrix(doc_term) freq_matrix <- data.frame(ST = colnames(temp_frequency), Freq = colSums(temp_frequency)) rm(temp_frequency)

# Convert to json

json_lda <- LDAvis::createJSON(phi = phi, theta = theta, vocab = vocab, doc.length = doc_length, term.frequency = freq_matrix$Freq)

return(json_lda) }

lda_json <- topicmodels_json_ldavis(fitted = desc_lda, corpus = myCorpus, doc_term = dtm)

# This chink does a multidimentional scaling intertopic distance map.
    # Required packages
    library(topicmodels)
    library(dplyr)
    library(stringi)
    library(tm)
    library(LDAvis)
    library(devtools)

Loading required package: usethis

devtools::install_github("cpsievert/LDAvisData")

Skipping install of 'LDAvisData' from a github remote, the SHA1 (43dd263a) has not changed since last install.
  Use `force = TRUE` to force installation

devtools::install_github("cpsievert/LDAvis")

Skipping install of 'LDAvis' from a github remote, the SHA1 (5067f7b8) has not changed since last install.
  Use `force = TRUE` to force installation

#' Transform Model Output for Use with the LDAvis Package
#' https://gist.github.com/trinker/477d7ae65ff6ca73cace
#' Convert a \pkg{topicmodels} output into the JSON form required by the \pkg{LDAvis} package.
#'
#' @param model A \code{\link[]{topicmodel}} object.
#' @param \ldots Currently ignored.
#' @seealso \code{\link[LDAvis]{createJSON}}
#' @export
#' @examples
#' \dontrun{
#' data("AssociatedPress", package = "topicmodels")
#' model <- LDA(AssociatedPress[1:20,], control = list(alpha = 0.1), k = 3)
#' LDAvis::serVis(topicmodels2LDAvis(model))
#' }
topicmodels2LDAvis <- function(x, ...){
    post <- topicmodels::posterior(x)
    if (ncol(post[["topics"]]) < 3) stop("The model must contain > 2 topics")
    mat <- x@wordassignments
    LDAvis::createJSON(
        phi = post[["terms"]], 
        theta = post[["topics"]],
        vocab = colnames(post[["terms"]]),
        doc.length = slam::row_sums(mat, na.rm = TRUE),
        term.frequency = slam::col_sums(mat, na.rm = TRUE)
    )
}
lda_json <- topicmodels2LDAvis(desc_lda)
serVis(lda_json, reorder.topic = TRUE)

Loading required namespace: servr
Failed with error:  ‘there is no package called ‘servr’’
If the visualization doesn't render, install the servr package
and re-run serVis: 
 install.packages('servr') 
Alternatively, you could configure your default browser to allow
access to local files as some browsers block this by default

library(RJSONIO)
x <- fromJSON(lda_json)
 

### this will map the tpic numbers from the rest of the analysis to the multidimentional scaling topic numbers.
print("Original - MDS#")

[1] "Original - MDS#"

for( i in 1:NumTopics) {
  print(paste(i, which(x$topic.order == i)))
}

[1] "1 16"
[1] "2 17"
[1] "3 12"
[1] "4 13"
[1] "5 18"
[1] "6 15"
[1] "7 7"
[1] "8 14"
[1] "9 4"
[1] "10 11"
[1] "11 8"
[1] "12 5"
[1] "13 9"
[1] "14 1"
[1] "15 2"
[1] "16 6"
[1] "17 3"
[1] "18 10"

# Generate bigrams

# Bigram setup
library(tidyr)
#nsf_funding_bigram %>%
#  count(bigram, sort = TRUE)
# remove bigrams with stopwords

bigrams_separated <- nsf_funding_bigram %>%
  separate(bigram, c("word1", "word2"), sep = " ")

bigrams_filtered <- bigrams_separated %>%
  filter(!word1 %in% stop_words$word) %>%
  filter(!word2 %in% stop_words$word)#replace all of the original words in the nsf_funding with lemmas
# This also needs to be done for the bigrams where the substatution happens for each of the two columns in the bigram. ???
lemma_bigram_new <- lemma_unique
#replace word 1 with lemma
lemma_bigram_new <- left_join(bigrams_filtered,lemma_bigram_new, by = c("word1" = "word")) %>%
  mutate(word1 = ifelse(is.na(lemma), word_clean, lemma)) # %>%
# drop last three columns id, word1 and word2
lemma_bigram_new <- subset(lemma_bigram_new, select = c(id, word1, word2))
lemma_bigram_new <- left_join(lemma_bigram_new, lemma_unique, by = c("word2" = "word")) %>%
  mutate(word2 = ifelse(is.na(lemma), word_clean, lemma)) # %>%

lemma_bigram_new <- subset(lemma_bigram_new, select = c(id, word1, word2))
#nsf_funding <- subset()

# new bigram counts:
bigram_counts <- lemma_bigram_new %>% 
  count(word1, word2, sort = TRUE)

bigram_counts

NA

ala bernhardlearns.com To speed things up, we will not work on every word instance but filter out only the unique ones. If needed the results for the unique word instance can be easily mapped back to the larger data-frame. Some of the techniques I am going to present are quite computational expansive, so if you have a much larger data set, then maybe they are not feasible. I have added system.time() to the particular work-steps so that you can see and decide for yourself.

Further, we apply some basic cleaning:

removing the possessive ending: ’s
removing all words containing non alphabetic characters (depending on the task at hand this might be a bad idea - e.g., in social media emoticons can be very informative)
removing missing values

bigrams_united <- lemma_bigram_new %>%
  unite(bigram, word1, word2, sep = " ")
bigrams_united %>%
  count(bigram, sort = TRUE)

#

https://www.tidytextmining.com/ngrams.html

library(igraph)

bigram_graph <- bigram_counts %>%
  filter(n > 6) %>%
  graph_from_data_frame()

library(ggraph)
set.seed(2017)
ggraph(bigram_graph, layout = "fr") +
  geom_edge_link() +
  geom_node_point() +
  geom_node_text(aes(label = name), vjust = 1, hjust = 1)

#topicFileName <- gsub(".csv","bigraph.png", filename)
#png(topicFileName, width = 12, height = 8, units = 'in', res = 300)
png(paste(topicDir,"/bigraph.png", sep=""), width = 12, height = 8, units = 'in', res = 300)

ggraph(bigram_graph, layout = "fr") +
  geom_edge_link() +
  geom_node_point() +
  geom_node_text(aes(label = name), vjust = 1, hjust = 1)
dev.off()

quartz_off_screen 
                2 

TF-IDF Analysis of bigrams

bigram_tf_idf <- bigrams_united %>%
  count(id, bigram) %>%
  bind_tf_idf(bigram, id, n) %>%
  arrange(desc(tf_idf))

bigram_tf_idf

# Graph the bigrams. THIS WILL TAKE SOME TIME
# categories are compressed
p <- bigram_tf_idf %>%
  arrange(desc(tf_idf)) %>%
  mutate(bigram = factor(bigram, levels = rev(unique(bigram)))) %>% 
  group_by(id) %>% 
  top_n(5) %>% 
  ungroup %>%
  ggplot(aes(bigram, tf_idf, fill = id)) +
  geom_col(show.legend = FALSE) +
  labs(x = NULL, y = "tf-idf") +
  facet_wrap(~id, ncol = 2, scales = "free") +
  coord_flip()

Selecting by tf_idf

#topicFileName <- gsub(".csv","BigramTFIDFHistogram.png", filename)
#ggsave(topicFileName, height = 14, device = "png")
p
ggsave(paste(topicDir,"/BigramTFIDFHistogram.png", sep=""), height = 30, device = "png")

Saving 20 x 30 in image

#png(paste(topicDir,"/TermsLDAHistogram.png", sep=""), width = 1280, height = 1600)
while (!is.null(dev.list()))  dev.off()

#Top documents per topic
library(ggplot2)
library(ggrepel)

number_of_documents = 5 #number of top docs to view
title <- paste("LDA Top Documents for", NumTopics, "Topics")

top_documents <- lda_gamma %>%
  group_by(topic) %>%
  arrange(topic, desc(gamma)) %>%
  slice(seq_len(number_of_documents)) %>%
  arrange(topic, gamma) %>%
  mutate(row = row_number()) %>%
  ungroup() #%>%
  #re-label topics
  #mutate(topic = paste("Topic", topic, sep = " "))

#topicFileName <- gsub(".csv","top_documents.csv", filename)
write.csv(top_documents, paste(topicDir,"/top_documents.csv", sep=""))

top_bigrams_for_docs <- bigram_tf_idf %>%
  group_by(id) %>%
  arrange(id, desc(tf_idf)) %>%
  slice(seq_len(number_of_documents)) %>%
  arrange(id, tf_idf) %>%
  ungroup()  %>%
  select(id, bigram) 
 
# make the 5 bigrams into 5 columns
top_bigrams_for_docs_wide <- 
  dcast(setDT(top_bigrams_for_docs), id~rowid(id, prefix="bigram"), value.var="bigram")   

# add the bigram columns to top_docs
### error - repeats of bigrams for too many documents....
# maybe not an error... documents in the smae topic are getting many of the same bigrams.
top_documents <- left_join(top_documents, top_bigrams_for_docs_wide,
            select(bigram1, bigram2, bigram3, bigram4, bigram5),
            by = c("document" = "id"))

#Join funding amount and title to top_docs_titles

nsf_funding_titles <- data_frame(document = MyData$AwardNumber,
                          title = MyData$Title, amount = MyData$AwardedAmountToDate)
# Convert the document #/grant number to a charctaer sting  
nsf_funding_titles$document <- as.character(nsf_funding_titles$document)

top_documents <- full_join(top_documents, nsf_funding_titles, by = "document")
#top_docs_titles <- merge(top_documents, nsf_funding_titles, by="document", all.x = TRUE)
#top_docs_titles <- subset(top_docs_titles, select=-c(row))

#top_documents <- full_join(top_documents, top_docs_titles, by = "document")
# add topic worth and rank
top_documents <- full_join(top_documents, topic_worth, by = "topic")

#topicFileName <- gsub(".csv","top_documents.csv", filename)
write.csv(top_documents, paste(topicDir,"/top_documents.csv", sep=""))

#top_docs_titles

#Top documents per topic


#FIGURE IS BAD AND NEEDS WORK, NOT REALLY A USEFUL VISUALIZATION ANYWAY
title <- paste("LDA Top Documents for", NumTopics, "Topics")

#define some colors to use throughout
my_colors <- c("#E69F00", "#56B4E9", "#009E73", "#CC79A7", "#D55E00", "#D65E00")

#customize ggplot2's default theme settings
#this tutorial doesn't actually pass any parameters, but you may use it again in future tutorials so it's nice to have the options
theme_lyrics <- function(aticks = element_blank(),
                         pgminor = element_blank(),
                         lt = element_blank(),
                         lp = "none")
{
  theme(plot.title = element_text(hjust = 0.5), #center the title
        axis.ticks = aticks, #set axis ticks to on or off
        panel.grid.minor = pgminor, #turn on or off the minor grid lines
        legend.title = lt, #turn on or off the legend title
        legend.position = lp) #turn on or off the legend
}

#customize the text tables for consistency using HTML formatting
my_kable_styling <- function(dat, caption) {
  kable(dat, "html", escape = FALSE, caption = caption) %>%
  kable_styling(bootstrap_options = c("striped", "condensed", "bordered"),
                full_width = FALSE)
}

word_chart <- function(data, input, title) {
  data %>%
  #set y = 1 to just plot one variable and use word as the label
  ggplot(aes(as.factor(row), 1, label = input, fill = factor(topic) )) +
  #you want the words, not the points
  geom_point(color = "transparent") +
  #make sure the labels don't overlap
  geom_label_repel(nudge_x = .2,  
                   direction = "y",
                   box.padding = 0.1,
                   segment.color = "transparent",
                   size = 3) +
  facet_grid(~topic) +
  theme_lyrics() +
  theme(axis.text.y = element_blank(), axis.text.x = element_blank(),
        #axis.title.x = element_text(size = 9),
        panel.grid = element_blank(), panel.background = element_blank(),
        panel.border = element_rect("lightgray", fill = NA),
        strip.text.x = element_text(size = 9)) +
  labs(x = NULL, y = NULL, title = title) +
    #xlab(NULL) + ylab(NULL) +
  #ggtitle(title) +
  coord_flip()

}

facet_wrap(~ topic, ncol = 10, scales = "free")

<ggproto object: Class FacetWrap, Facet, gg>
    compute_layout: function
    draw_back: function
    draw_front: function
    draw_labels: function
    draw_panels: function
    finish_data: function
    init_scales: function
    map_data: function
    params: list
    setup_data: function
    setup_params: function
    shrink: TRUE
    train_scales: function
    vars: function
    super:  <ggproto object: Class FacetWrap, Facet, gg>

word_chart(top_documents, top_documents$document, title)

NA
NA

#Distribution of probabilities for all topics

ggplot(lda_gamma, aes(gamma)) +
  geom_histogram() +
  scale_y_log10() +
  labs(title = "Distribution of probabilities for all topics",
       y = "Number of documents", x = expression(gamma))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Transformation introduced infinite values in continuous y-axis
Warning: Removed 5 rows containing missing values (geom_bar).

#Probability distribution in each topic

ggplot(lda_gamma, aes(gamma, fill = as.factor(topic))) +
  geom_histogram(show.legend = FALSE) +
  facet_wrap(~ topic, ncol = 4) +
  scale_y_log10() +
  labs(title = "Distribution of probability for each topic",
       y = "Number of documents", x = expression(gamma))

Error in ggplot(lda_gamma, aes(gamma, fill = as.factor(topic))) : 
  could not find function "ggplot"

https://knowledger.rbind.io/post/topic-modeling-using-r/ We can correlate the topic by the metadata available, here will will use the predefined category of the documents. I construct a dataframe binding the Document Number, theta for the document, and its associated category. The column means of theta, grouped by category is calculated.

theta <- posterior(desc_lda)$topics %>% as.matrix
x <- as.data.frame(row.names(theta), stringsAsFactors = FALSE)
colnames(x) <- c("GrantNumber")
x$LessonId <- as.numeric(x$LessonId)

Error in `$<-.data.frame`(`*tmp*`, LessonId, value = numeric(0)) : 
  replacement has 0 rows, data has 38

Financial analysis on Award Amount

#Join abstracts to funding amount
library(dplyr)

nsf_amount <- data_frame(id = as.character(MyData$AwardNumber),
                            title = MyData$AwardedAmountToDate)
nsf_amount

#Continue join

nsf_amount <- nsf_amount %>% 
  unnest_tokens(word, title) %>% 
  anti_join(stop_words)

Joining, by = "word"

desc_tf_idf <- full_join(desc_tf_idf, nsf_amount, by = "id")
desc_tf_idf

# 
TopicInvestmentBigram <- full_join(top_documents, topic_worth, by = "topic")
write.csv(TopicInvestmentBigram, paste(topicDir,"/TopicInvestmentBigram.csv", sep=""))

Identify the top element codes for a topic. Each grant has a gamma associated with each topic. Each grant has one or more element code. Program codes indicate cross-NSF initiatives. The element code is home but funding also arrises from the cross NSF programs. For a first aproximation analaysis calculate the element contribution to a topic by multiplying the grant amount for each grant in the for each topic ProgramInvestment = Sum of gamma * amount divided evenly among elment codes. This distribution is an approximation because the split for funding may not be 50/50. eg. if there is one element code all spending is assigned to that one element. If there are two the investement will be split between them. Money may actually be coming form program codes as well.. Output Data Structure: Topic, Element Code, Total Investment